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Abstract 

By directly integrating the Schrodinger equation starting in the transmission region and working 
backwards through the barrier, the tunneling probability can be determined for arbitrary potential 
barriers. The method employs techniques familiar to undergraduates and is used here to study 
resonant tunneling. 
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I. INTRODUCTION 



Quantum tunneling is a favorite conceptual topic for students. It is a notion of something 
so very different from what is expected classically that can be described so easily by invoking 
memories of throwing balls at walls. Students are encouraged to find connections with 
frustrated total internal reflection in order to further cement their understanding of matter 
as waves. Both in optics and quantum mechanics instructors can note how the continuity of 
the wave function and its derivative across boundaries implies that the wave cannot abruptly 
go to zero. This enables students to see why a (typically small) portion of the wave can 
tunnel through a barrier. 

On the other hand, the quantitative aspects of tunneling are a different story for students. 
As usual, very simple situations can be calculated by hand like the square barrier but more 
typical barriers found in the lab require the use of a computer and an algorithm that can 
apply the conceptual physics they have learned (the boundary conditions) in an iterative 
fashion. 

In this work present a different algorithm to calculate the tunneling probability of a 
particle with known energy through an arbitrary one- dimensional potential barrier. It is both 
fast and accurate but it also uses tools that most students are familiar with in their studies 
of the Schrodinger equation. Specifically it involves the direct integration of the Schrodinger 
equation in a manner very similar to the shooting method employed to find the eigenstates 
of an arbitrary potential well. However, instead of needing to adjust parameters to find a 
particular eigenstate, students can directly inspect the results for any given particle energy 
and determine both the tunneling probability and the shape and nature of the wavefunction 
inside the barrier. 



II. OTHER METHODS 

The most common approach to calculating tunneling probabilities is to consider the 
barrier to be a collection of square barriers. In the WKB approach, only the exponen- 
tially decaying portion of the wavefunction is kept and integrated through all the slices.- 
In the matrix transfer method, the boundary conditions among all the slices are carefully 
calculated." Specifically, at every boundary between the square slices the wavefunction 

2 



and its slope are continuous. In each slice the wavefunction is composed of two components: 
either a right and left traveling wave with a wavelength determined from the kinetic energy 
(the difference between the total energy and the barrier height); or a growing and decaying 
exponential whose growth rate is determined from the (negative) kinetic energy. Often these 
boundary condition equations are described in a matrix formalism as they are simple linear 
equations relating the incoming and outgoing wavefunctions along with the barrier heights 
of the slices. The effect on the incoming wave by the barrier is then modeled by a single 
matrix that can used to solve for the tunneling probability. 

There are also some approaches in the literature that have more directly integrated the 
Schrodinger equation but all do a forward propagation as opposed to the backward one 
described below .-^ These approaches use both numeric and analytical methods to determine 
the phase of the incoming wave that enables solely a right-traveling wave in the transmission 
region. The method below does not require such adjustments and simply gives both the 
wavefunction in the tunneling region and the tunneling probability after a single direct 
integration. 



III. METHOD 

Consider a tunneling situation as laid out in Figure [TJ The first and third regions have a 
constant potential while the middle region can have any form, including discontinuities and 
regions where the particle is classically allowed. Region I can have right- and left-traveling 

waves 

^ = Ae iklX + Be~ iklX (1) 
while Region III only has a right traveling wave 

tjj ni = Fe ikmX (2) 

where 

h = sJ{2m 2 /h 2 )E (3) 

and 

km = ^(2myh 2 ) (E-V m ). (4) 

Using a fourth-order Runge-Kutta technique we numerically integrate the real and imag- 
inary parts of the Schrodinger equation from the right edge of Region II (x = L) to the left 
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edge {x = 0). Note that since the Schrodinger equation does not have any single deriva- 
tives in it a Numeroff approach can also be used.- Note also that in Mathematica you can 
integrate complex numbers with just one call to the Runge-Kutta solver (NDSolve). Since 
both the wavefunction and its slope will be the same on both sides of the boundary between 
Regions II and III, the initial conditions are determined by arbitrarily setting F = 1 and 
using the form from Region III: 

4>{L) = e ikmL and ij)'{L) = ik m e ikluL . (5) 

To determine the transmission probability, T, we need to find the value of A: 



T = — 
h 



hU 1 (6) 



h \a\ 2 ' 

This is done by investigating the value of the wavefunction and its slope at x = where, 
according to Eq. [TJ 

ip(0) = A + B and ^'(O) = iki(A - B). (7) 

Once again we have used the equality of the value and slope of the wavefunction across a 
boundary. 

Combining Eqs. |6] and [7] yields 

T _ fc ni 4 

h |^(0) -#'(0)Ai| 2 ' 

Once the Schrodinger equation has been numerically integrated, the transmission probability 
is easily calculated. 

This method employs many techniques used when teaching the numerical solution of 
eigenstates for arbitrary energy wells. In those situations students are taught to employ the 
shooting method to find energies that produce physically allowable wavefunctions. Students 
start with an initial value and slope (these can usually be pretty arbitrary), typically on one 
side of the energy well to be considered, and use a numeric integrator to predict the shape 
of the wavefunction in the well. They then adjust the energy involved, which changes the 
shape, and they reject any energies that cause the wavefunction to asymptotically approach 
infinity at the other side of the well. In this new application the initial value and slope are 
set by Eq. |5]and the solution at every energy is used to calculate the transmission probability 
using Eq. |HJ The two major differences in this new application are that the integration is 



done backwards spatially and that both the real and imaginary parts need to be integrated, 
as is illustrated in Figure [2j If you only do the real part (as is often done in the shooting 
method application) you are unable to calculate the transmission coefficient as seen in Eq. |HJ 



IV. EXAMPLES AND COMPARISONS 

The transmission coefficient (T) as a function of particle energy for the potential shown 
in Figure [1] is given in Figure |3j The top curve is the result of the current method while 
the lower lines use the transfer matrix method with varying number of slices of Region II. 
Ultimately both approaches converge to the same result at every energy. 

It is interesting to compare the transfer matrix method with the new method where 
the number of slices is compared with the number of steps that the Runge-Kutta method 
employs. The transmission probability versus energy for the arbitrary barrier shown in 
Figure [1] is given for the total step number ranging from 5 to 12 in the inset of Figure [31 
The curve with 300 steps is also shown. It is clear that the number of steps needed for 
the Runge-Kutta method is far less than the number of slices needed in the transfer matrix 
method to achieve the same accuracy. Note, however, that one should really compare the 
number of calculations involved when doing these comparisons. A fairer comparison would 
need to multiply the number of Runge-Kutta steps by four, though this still shows that the 
current approach compares favorably to the tranfer matrix method. 

V. RESONANT TUNNELING 

As an example to show the pedagogical uses of the current method, we consider resonant 
tunneling. Specifically we compare the wavefunction in the barrier region to the eigenstates 
expected for a simply-shaped barrier. 

Consider the potential shown in the inset of Figure HI This parabolic potential barrier is 
given in electron volts by 



10(x- 1) 



2 if < x < 2 



(9) 



V(x) 







otherwise 
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The analogous potential well that is not truncated has resonant energies at 

E n = (n + = (n + ~j 1.232 eV where n = 0, 1, 2, . . . . (10) 

The transmission probability as a function of energy is shown in Figure HJ The resonance 
peaks shown correspond very nearly with the eigenenergies of a parabolic well. At lower ener- 
gies where the resonance peaks are very sharp the energies are the same as the eigenenergies. 
As the peaks become broader, the resonant energies become larger than the eigenenergies 
by as much as 20% for n = 11. 

The reason the resonant energies grow larger than the eigenenergies as the energy in- 
creases is due to the boundary conditions that the wavefunction has to match at x = 
and x = 2. This can be seen by comparing the resulting wavefunctions (both the tunneling 
wavefunction and the eigenfunction for the parabola) as seen in Figure El Close inspection 
of the wavefunctions near the boundaries indicates the differences between the tunneling 
wavefunction and the eigenstate. While the eigenstate is decaying to zero in all cases, the 
tunneling wavefunction is forced to match the boundary condition at the right edge. At low 
energies there is little difference as the exponential rise is very steep but at higher energies 
there is a higher bend needed that explains the rise in the energy compared to the analogous 
eigenenergy. 



CONCLUSION 



We have discussed a new method for calculating both transmission probabilities and 
wavefunctions for a particle tunneling through an arbitrary one-dimensional barrier. The 
approach is applicable at the undergraduate level as it uses common tools related to the 
shooting method for finding potential well eigenstates. It is fast and accurate and enables 
the study of complex phenomena like resonant tunneling. 
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FIG. 1. Arbitrary potential barrier. The electron approaches in Region I and tunnels through 
Region II into Region III. 
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FIG. 2. The real (solid) and imaginary (dashed) parts of the tunneling wavefunction for the barrier 
shown in Figure Q] for a particle energy of 7.5 eV. This yields a tunneling probability of 93%. 
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FIG. 3. The transmission probability, T, as a function of particle energy for the potential barrier 
shown in Figure [TJ This compares the current method with the matrix method with the number 
of slices set to 10, 20, 50, 100, 200, 500, 1000. The current method is the top curve. In the inset 
the current method is limited to a total number of steps ranging from five through twelve. 
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FIG. 5. Comparison of the resonant tunneling wavefunction with the eigenfunction for a parabola 
for n = 5 (top) and n = 9 (bottom). The probability density (ijj*ip) is plotted as a function of 
position inside the potential barrier (shown in Figure 0]). The curves are all normalized to a max 
of unity. The resonant wavefunction has higher energy so that the curve can bend to better meet 
the boundary conditions. 
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